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We examine the stability of magnetic order in a classical Heisenberg model with quenched random 
exchange couplings. This system represents the spin degrees of freedom in high-Tc compounds with 
immobile dopants. Starting from a replica representation of the nonlinear cr-model, we perform 
a renormalization-group analysis. The importance of cumulants of the disorder distribution to 
arbitrarily high orders necessitates a functional renormalization scheme. From the renormalization 
flow equations we determine the magnetic correlation length numerically as a function of the impurity 
concentration and of temperature. From our analysis follows that two-dimensional layers can be 
magnetically ordered for arbitrarily strong but sufBciently diluted defects. We further consider the 
dimensional crossover in a stack of weakly coupled layers. The resulting phase diagram is compared 
with experimental data for La2-a;Sra:Cu04. 



PACS numbers: 75.10.Nr, 74.72.-h, 75.50.Ee 



I. INTRODUCTION 

Although the interest in disordered spin systems 
reaches back several decades (for a review, see e.g. Ref. 
|l|), interest has strongly revived in recent years by the 
recognition that high-Tc compounds exhibit phases with 
antiferromagnetic (AFM) order, spin-glass order or stripe 
order in certain ranges of temperature and doping. A 
prominent example for such materials is La2-a:Sra;Cu04, 
in which every Cu atom carries a spin i (for an overview 
on this material, see e.g. Ref. ||). 

In these materials, the spins are located in weakly 
coupled layers. Within each layer, the spins constitute 
a square lattice with antiferromagnetic exchange cou- 
phng between nearest neighbors In the undoped system, 
these spins can be represented to a good approximation 
by a classical two-dimensional model since quantum fluc- 
tuations lead to a merely auantitative renormalization 
of the classical parameterso Although a weak coupling 
between the layers and a weak easy-plane anisotropy 
are present in these materials (both are about five or- 
ders of niagnitude smaller than the isotropic intra-plane 
exchanged) , they become relevant only on relatively large 
scales. Thus, on the finite length scales of experimental 
relevance, the spins can be described in a first approxi- 
mation by a classical Heisenberg antiferromagnet in two 
dimensions. 

Doping induces holes in the layers, which can lead to 
quenched frustration of the exchange couplings. FrustEtUj 
tion can occur when the holes are localized individuallyQ'El 
as well as when they condense into a topologically defec- 
tive array of stripeS|-»diich act as antiphase domain walls 
for the AFM order.laa Both cases can be represented by 
an effective bond-disordered Heisenberg model. Such a 
system is expected to display generic spin-glass behavicp, 
as is observed in experiments beyond a critical doping.t2l 

On the theoretical side, Heisenberg spin glasses are 
much less understood than Ising spin glasses. While it 
is well established for the pure two-dimensional system 
that the magnetic correlation length ^ decays exponen- 



tially with increasing temperature,lli the dependence of 
^ on disorder is controversial. At T = it is not clear 
whether ^ is finilfi-foiipiery small but finite concentration 
of defect bondsBll3lijO or whether ^ is infimte up to a 
critical concentration even for strong defectsJl3 Concern- 
ing the temperature dependence-of ^, thep, is no consen- 
sus as to whether it is reentrantli3 or not. Ell To readdress 
these unsettled issues, we develop a theoretical approach 
for bond disordered spin glasses which complements pre- 
vious approaches. In particular, we calculate the corre- 
lation length and compare our results to the results of 
previous theoretical approaches and to experiments on 
high-Tc compounds. 

The outline of the paper is as follows. In Sec. II we 
establish the model and briefiy review the mechanism 
by which holes generate magnetic textures that provide 
the nucleus for a reduction of magnetic order. Based 
on this mechanism, we further motivate why it is de- 
sirable to develop an alternative to previous theoretical 
approaches. In Sec. HI, we derive the functional renor- 
malization group equations. These fiow equations are 
evaluated numerically and our results for the magnetic 
correlation length are compared to experimental data in 
Sec. IV , where we also consider the dimensional crossover 
from 2D to 3D behavior. We conclude with a discussion 
of our results in Sec. 0. 



II. MODEL 

In a system with sparse substitutional doping, the 
dopants can be considered as randomly distributed with 
negligible correlations and as quenched over a wide tem- 
perature range. In doped Mott insulators, the effect of 
the dopants is twofold: they induce holes in the cuprate 
planes and at the same time they provide a random po- 
tential that localizes these holes. Therefore, they are 
commonly assumed to be localized iQgrtjiie oxygen atoms 
between neighboring copper atoms.Ql3ll3 

A simple argument suggests that such a hole trans- 
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mutes the antiferromagnetic superexchange coupling be- 
tween the neighboring copper spins into a ferromagnet: 
one. Thereby frustration is induced in the spin system 
The irrelevance quantum fluctuations of the spins in 
the pure systemQ (in the sense that they lead only to a 
quantitative renormalization of classical parameters such 
as the spin stiffness) suggests that they may be Jifiglected 
in a good approximation also for small dopingJHI 



appendix The ground state for a single defect with 
A > Asingio remains planar for all N > 2, i.e., apart from 
a global rotation, the defect texture for Heisenberg spins 
{N — 3) is identical to the texture for XY spins {N — 2). 
Far away from the defect bond, the texture is described 
by the solution of the Laplace equatiorO 



(4) 



A. Definition 

Starting from a classical description of the AFM, we 
may examine an equivalent model where all spins of one 
bipartite sublattice are flippied and the sign of all ex- 
change couplings is reversed,c3 i.e., we now consider the 
coupling of the undoped system to be ferromagnetic and 
the defect couplings to be antiferromagnetic. We base 
our analysis on the Hamiltonian 



(1) 



for spins on a (hyper)cubic lattice in d dimensions (a 
single layer is described by d = 2). Spins are treated 
as classical A^-component vectors of unit length, (r) — 
1. A spin at site r is coupled to its nearest neighbor 
in direction i = 1, ■ ■ ■ ,d via the exchange coupling Ji(r). 
We define ViS{r) := S{r+ai) — S{r) for a basis vector a^. 
The global symmetry 0{N) of spin rotations is preserved 
for arbitrarily disordered exchange couplings. 

It is convenient to rewrite the exchange coupling as 



J,{r) = [1- A,(r)]J 



(2) 



with the value J > of the pure system and the quenched 
random variable Aj;(r). Frustration effects emerge for 
Ai(r) > 1, when a ferromagnetic bond becomes antifer- 
romagnetic. To further specify the nature of disorder, we 
assume a bimodal distribution of the exchange couplings. 



A,:(r) 



A with probability p, 
with probability 1 



(3) 



For square lattices as in La2-xSra;Cu04, the concentra- 
tion X is related to this probability by a; = 2p. Correla- 
tions between different bonds are assumed to be absent. 



B. Single defects 

Before we address the magnetic order in the presence 
of a finite concentration of defects, it is instructive to 
recall briefiy the physics of a single defect for A > 2. 
While the ground state of the pure system is coUinear 
(all spins are strictly parallel), it is canted (i.e., no longer 
coUinear) for a r^Pftl©] defect beyond a certain critical 
strength Asingic-"ilj'c3 The threshold Asingi(r-T= d value 
can be obtained from a spin-wave calculationEil (see also 



and the spins approach a coUinear configuration, 5(r) — » 
So for r —>■ oo. To be specific, we discuss a defect bond 
located at r = and oriented in direction a. It acts as a 
source of a dipolar distortion with a moment /2 perpen- 
dicular to So- For d = 2, in particular. 



S{r) - So 



fi a - r 



(5) 



The amplitude of the dipole moment is determined by the 
nonlinearity of the model, which may be considered as in- 
teraction between spin waves. The magnetic textures of 
several defect bonckpare subject to a dipolar interaction 
at large distances.E3o The effects of a finite density of 
defects are nontrivial because of the frustrated long-range 
interaction between the dipoles. In the presence of more 
than two defects, the interaction among the dipoles typ- 
ically becomes frustrated. Then the ground state is no 
longer planar for A > 2. 



C. Finite defect concentration 

Depending on the relative position of several defect 
bonds in a cluster, the generated texture may have dipo- 
lar or higher-order multipolar structure. In contrast to 
electrodynamics, the multipolar moments of defect tex- 
tures are not additive since the spins act as a nonlinear 
medium for the dipoles, giving rise to many-body interac- 
tions. In particular, clusters of defects can lead to cant- 
ing already for defect-bond strengths below Agingic, i-e-, 
in clusters of defects the threshold strength of the defect 
bonds is reduced. The threshold values of some specific 
bond configurations in the XY model are given in Refs. 
pl| and |2^ (cf. in particular Table I in the latter one). 
In general, threshold values in the XY system are upper 
bounds for the thresholds in a Heisenberg system since 
with increasing A the spins have a larger space of canted 
states which they may explore to minimize energy. 

Thus, as soon as the defects are antiferromagnetic 
(A > 1) they can induce canted teidiips-due to the 
presence of arbitrarily large clustersEjOSE^ Hence the 
ground state will be canted for every p > and A > 1 
(assuming that the defect distribution has no patholog- 
ical spatial correlations). Therefore one may expect a 
system with a certain density of weak defects 1 < A < 
Asingic to be qualitatively equivalent to a system with a 
(possibly much) smaller concentration of strong defects 

A > A.inglo. 
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The presence of canting certainly implies a reduction 
of magnetization. However, one cannot draw direct con- 
clusions about the range of magnetic order. For a suffi- 
ciently high density of strong defects, one certainly has 
short-range order. For defects of a weaker strength or a 
lower concentration, one can have also quasi-long-range 
order (p^riJl a XY model with quenched and uncorrelated 
dipolegE3o). A phase with true long-range order is also 
possible but unlikely since it would require a highly or- 
dered dipole configuration. I_. 

It is instructive to recall a "duality" relatione^ in the 
{p, A) parameter space, since in a mixture of ferromag- 
netic and antiferromagnetic bonds there are two ways to 
declare one type as "regular" bond and the other type as 
"defect" bond (see appendix 0) . This duality implies a 
relation 



lAFM 



Afm(1 -p) 

{l-p)-l 



(6) 



between boundary lines Afm(p) and Aafm(p) limiting 
regions with ferromagnetic order (for A < Afm) and 
antiferromagnetic order (for A > Aafm) in the (p, A) 
plane. Since the duality relation maps the region pA > 1 
onto pA < 1 it is sufficient to examine the latter region. 

To perform first qualitative estimates, one may con- 
sider the disorder-averaged exchange coupling, which is 



Mr) = (l-pA)J 



(7) 



for the bimodal distribution (g). Thus, there is a ten- 
dency towards the formation of ferromagnetic order for 
pA < 1 (where J > 0) and a tendency towards antiferro- 
magnetic order for pA > 1 (where J < 0). Certainly, the 
presence of order requires more than such a tendency. A 
further minimum requirement should be that the relative 
fluctuations 



J2(r)-J,(r) _p(l-p)A2 



Mr) 



(l-pA)2 



(8) 



of the exchange couplings must be small. Since A > 1 
and pA < 1 in the range of interest, ^ 1 requires 
pA^ < 1. A very crude estimate of the boundary from 
cr-1 suggests 



Afm{p)-P''^^- 



(9) 



For the special case A = 1, the defect bonds have 
a vanishing exchange coupling and the system is bond- 
diluted. The presence of magnetic order (in the sense of 
a finite magnetization) then requires that a finite fraction 
of spins is connected by regular bonds. This is the ca^ 



below the percolation transition, i.e. for p < ^ in = 2.1 



Therefore one expects Afm (p) — 1 for p > 



D. Previous work 

There are only a few approaches in the literature aim- 
ing at a more sophisticated analysis, which we briefly 



summarize in order to highlight our motivation to re- 
consider this problem in an alternative way. First, the 
coherent-potential approximation (CPA) provides a sim- 
ple self-consistent approach to determine an effective spin 
stiffness for the random system. For N = 2, the CPA 
yieldscZl a transition line Ap^'^ (p) smoothly interpolating 
between Ap^^(O) = d (reflecting the canting threshold 
for individual defect bonds) and Ap^{p) = 1 for p > i 
in d = 2 and for p > | in rf = 3. In d = 2, the location 
of the percolation transition is captured exactly. This 
has to be considered as a fortunate coincidence which 
is absent in dr-^^ 3, where the location is found only 
approximatelyED In the limit — > oo thCj-CPA yields 
that the stability of order requires < l,t3 which re- 
sults in Ap^{p) « p~^/^ for small p in agreement with 
the naive estimate (|9|). 

Unfortunately, jtlie interpretation of the CPA in some- 
what ambiguous.c3 According to its construction, the 
CPA replaces the disordered system by a ferromagnetic 
one with a homogeneous effective exchange constant. 
The strength of this effective coupling is determined from 
a selfconsistency condition that considers the fluctuations 
of a single bond. Thus, the CPA misses effects of clusters 
of defect bonds. The putative transition line is defined 
as the border line up to which selfconsistcnt solutions ex- 
ist. A priori, this border line can be interpreted in two 
ways: either as the onset of the canting instability or as 
the onset of short-range magnetic order. Therefore, the 
CPA can be only of qualitative use. 

To gain insight into-tiie nature of the ground state, 
Gawiec and GrempelH3'c3 (GG) performed numerical 
studies of the XY model in c? = 2. Their data suggest 
a transition line Ap^(p) between a phase with quasi- 
long-range order and short-range order. It starts from at 
the canting instability Ap^(O) — Agingie ~ 2 and follows 
Ap^(p) = 1 beyond the percolation transition. Thus, 
the data suggest short-range order for an arbitrarily small 
concentration of defects of a strength exceeding the cant- 
ing instability of single defects in agreement with CPA. 
However, as stated in Ref. |2^, the question of whether 
one can have order for A > Agingio at sufficiently small p 
remains far from settled due to finite-size effects. 

For the Heisenberg model with the bimodal bond dis- 
tribution, there are only very few numerical studies. For 
the special case of two, dimensions and A = Agjngic — 2, 
Nonomura and OzekiE^ postulate order for p < 0.11 from 
an exact-diagonalization method for S* = 1/2 quantum 
spins. This implies an even larger stability of the classi- 
cal model, and in particular Afm(O) > 2. However, also 
their conclusion has to considered with care because of 
finite-size effects. i— ii— i 

Previous analytic approacheai3lij are implicitly re- 
stricted to defect bond strengths exceeding the canting 
threshold, A :s> Agingio- The spin system with ran- 
dom bonds is replaced by a spin system with homoge- 
neous bonds coupled to an additional canting field that 
also generates dipolar spin textures. In order to shed 
some light on the quality of this replacement, we perform 
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a Hubbard-Stratonovich transformation introducing an 
auxiliary bivectorial field fi via 

eM-H[S]/T) = J cxp[-H[5,/]/T],(10a) 



with 



H[S,f] = jY,S^^[VJir)r - Mr) -VJir) 



(10b) 



In this representation, S and fi are thermally fluctu- 
ating variables, S with spheri cal c onstraint, fi uncon- 
strained. In the Hamiltonian (10b) of the transformed 
system, the spins interact directly via the homogeneous 
exchange coupling J. In addition, they couple also to the 
canting field fi. It is precisely such a coupling that was 
considered in Refs. [T^Jl^ . 

The Hubbard-Stratonovich transformation - which 
yields an exact representation of the original model (|l|) - 
shows that the field fi has a local selfinteraction potential 



yy(r,r') =(5y(5r^r'J/A,;(r). 



(11) 



On regular bonds with Ai(r) = 0, fi{r) is suppressed. On 
defect bonds with Ai(r) > 0, /i(r) has finite fluctuations. 
For ViS'(r) = one would have 



T 



(/r(r)/j'(r'))=<5^%,r'jA,(r). 



(12) 



However, it is crucial to retain the full correlations be- 
tween the fluctuations of S and fi. The spin- wave saddle- 
point equation 



be a reason for which GI find a reentrant temperature 
dependence of 

On the other hand, Cherepanov et alS3 consider fi as 
a quenched field with Gaussian correlations. This is in 
contradiction to the annealed nature of fi as is revealed 
by the transformed model (p^). In addition, spin rota- 
tion symmetry is explicitly broken. This treatment is 
based on the assumption that the spin textures freeze 
at low temperatures. However, a spontaneous breaking 
of this symmetry piti d < 2 is ruled out by the Mermin- 
Wagner theorem.E2l Thus, one may worry that disorder 
effects may be overestimated also in this approach. Pos- 
sibly, the artificial symmetry breaking is related to the 
spurious generation of random fields-, in a replica treat- 
ment as found by Cherepanov et alS-1 



III. RG ANALYSIS 

Although the Hubbard-Stratonovich transformation 
was useful to relate previous work to the original model, 
the introduction of / is accompanied with additional dif- 
ficulties. Therefore we choose to continue to work with 
the original model where the field / is integrated out. Be- 
fore we perform a renormalization-group (RG) analysis, 
we use the standard replica trick to treat disorder. 



A. Replica representation 



After n-fold replication the Hamiltonian reads 



V,5(r) ^ /,(r) 



(13) 



with the abbreviation 

n 

Mr) ■.= J2i^J''{r)? 



(14a) 



(14b) 



immediately shows that fi induces a canting of the spin 
field. If fi{r) is considered as fixed and nonvanishing only 
on a single bond, it induces a dipolar spin texture. It is 
important to notice that the original spin rotation sym- 
metry of Hamiltonian (p]) is preserved in the transformed 



model ( 10b ) only if fi is rotated simultaneously with S. 

Instead of solving the full problem of two fluctuating 
fields, Refs. ^ 14 proceed with additional assumptions 
about the nature of fi. Glazman and Ioselevicht3 (GI) 
consider a Hamiltonian, where fi (r) has a fixed length on 
the defect bonds. Only the orientation of fi{r) is consid- 
ered as a thermal degree of freedom. Thereby the rota- 
tion symmetry is preserved. However, fixing the length 
of fi(r) contradicts to Eq. ( p2[ ) which shows that also 
the magnitude of fi{r) is a thermally fluctuating quan- 
tity that vanishes in the limit of zero temperature. This 
means that in the approach of GI the strength of dis- 
order is overestimated at low temperatures. This may 



(Upper Greek indices label replicas.) Since we assume 
that the probability distribution of Ai(r) is uncorrelated 
and identical for all bonds, disorder averaging leads to the 
local and translation symmetric Hamiltonian {H includes 
the factor 1/T) 



n = Y,f^^KMr)~RiMr)) 



The cumulant function R is specified by 



R{Mr)) :=ln exp[A,(r)i^V^(r)/2] 



(15) 



(16) 



and depends implicitly on K 
case of bimodal disorder (^) 



i?(^) 



In 



J/T. For the special 
(17a) 



1 — p + j3 e 

yAK^ forAifV^O, , . 

iAi^V + lnp forAXV^oo. ^ ' 



5 



Note that naturally = in the absence of disorder 
(for p = or A = 0) and that R{ip) = ^AK^p for p = I, 
which amounts to an unfrustrated dual model with a cor- 
responding stiffness (1 — A) J. In the general case, R{ip) 
is a nonlinear function which has a linear asymptotics for 
small and large arguments, cf. Rq. (|l7|). 

Remarkably, the energy depends on ip only through 
the combination 



(18) 



A priori, it is not clear whether spin fluctuations are 
governed by the behavior of h(^jj) at large or at small 

The stability of the ferromagnetic state with respect 
to large-scale spin-wave deformations depends on h^ip) 
at small ip and requires h'{0) > 0, i.e., pA < 1. (For 
pA > 1 the dual antiferromagnetic ground state is locally 
stable.) This condition is equivalent to the requirement 
that the average exchange coupling (|^) should remain 
ferromagnetic. 

While the effective stiffness is defined from the slope 
of the function h{ip) at i/j — 0, its behavior h{ip) ~ 
(1 — A)Kil; for -0 ^ oo (which is equivalent to T ^ 0) 
reflects the presence of frustration. For A < 1 — i.e., 
in the absence of frustration — h{'ip) is positive for all ip 
and fluctuations can renormalize the stiffness to a smaller 
but positive value. For A > 1 — i.e., in the presence of 
frustration — h{tjj) is negative at large ip and fluctua- 
tions can renormalize the stiffness to zero, signaling the 
destruction of magnetic order. 

In the replica representation, the canting instability 
for a single defect can be retrieved easily. To this end, 
one has to keep the cumulant function ( [T7| ) only on a 
single bond and to switch it off on all other bonds. In 
the limit T — » (considering the statistical weight for 
fixed n and an fixed but arbitrary spin configuration), 
the replicas decouple since i?(V') « iAJsTV + Inp. The 
exchange coupling of the defect bond within each replica 
is (1 — A) J independent of p. Thus, canting occurs in the 
case p < 1 for A > Agingio as in the case p = 1. While 
the threshold is independent of p, the disorder-averaged 
dipole moment depends on p. Since ip becomes arbitrarily 
small slightly above the canting threshold, this dipole 
moment cannot be calculated from the large-?/' limit of 
h. 

Therefore it is in general important to retain the global 
fimctional form of h(ip) [or, equivalently, R{'ip)]. An ap- 
proximation of this function by a Taylor series near -ip ^ 
would amount to an expansion in cumulants of the de- 
fect distribution. The relevance of R{ip) for large ip shows 
that one could miss essential physics by dropping high- 
order cumulants and - in particular ~ by using a Gaus- 
sian distribution for A. Even worse: the function R{4') 
is nonanalytic, i.e., its behavior at large ^p is outside the 
radius of convergence of a cumulant expansion. 



B. Flow equations 

In order to address the question of magnetic order in 
the presence of general fluctuations oi 'ip, we generalize 
the renormalization-group analysis of Polyakovllll to the 
replicated model. 

Instead of working with the function hiip), it is more 
physical to use an effective bare stiffness and an effective 
cumulant function deflned by 



Ko := K -2R'{0) = {l-pA)K, 



(19a) 



RoW := RW - Vi?'(0) = RW - -pAK^{19h) 

as to remove the linear contribution of R, which repre- 
sents a trivial renormalization of the stiffness. The re- 
placement of the original quantities by the effective ones 
leaves the energy ( p^ ) invariant since — 2R{ip) = 
Koi> - 2i?o(^). „ 
Following Polyakov's analysis of the pure system,Eil 
we renormalize the model by a momentum-shell inte- 
gration. So far, the Hamiltonian (15) was written in 



a way explicitly retaining the lattice of the spin sites. 
We choose the lattice spacing as unit of length and go 
over to a continuum representation by smoothly interpo- 
lating the fleld S{r) between the lattice sites (preserv- 
ing the normalization S'^(r) = 1 everywhere). In this 
limit, partial spatial derivatives replace the differences in 
Mr)^J:l=l[^^S"{r)r and 



(20) 



The replacement of differences by partial derivatives and 
of the sum by an integral should be a reasonably good 
approximation, i.e., the replacement of ip should lead 
to relative errors small compared to unity for an arbi- 
trary spin conflguration (even with inhomogeneities on 
the scale of the lattice spacing). A small relative error 
in tp is equivalent to a small relative misrepresentation of 
J and/or A since these quantities enter the Hamiltonian 
only as a product with ip. We further use the approxi- 
mate replacement of the cubic Brillouin zone by a spher- 
ical one. Thus, to preserve the volume of the Brillouin 
zone, its radius A has to be fixed by A'' = d{2TT)'^/Sd 
with Sd := 27r'-'^/^^ /r{d/2) the surface of the d dimen- 
sional unit sphere. In d = 2 specifically, A^ = Att. In 
order to demonstrate that the fundamental frustration 
mechanism does not get lost due to the continuum ap- 
proximation, we briefly rederive the canting threshold in 
Appendix ^ 

According to the scheme of the momentum-shell renor- 
malization group, we now integrate out the spin modes 
with wave vectors in the shell Ae~'^^ < k < A. Thereby 
the original spin field S'" is mapped onto the slowly vary- 
ing background field s". They are related by 



5" = yr 



AaAa ^ 



(21) 
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with the vector fields forming a local orthonormal ba- 
sis {e^,-- - ,(^_i,s"} at each site r in each rephca a. 
We employ the sum convention for pairs of Latin indices 
(a,5=l,---,A^— 1) only. The field generates an in- 
finitesimal spin rotation and has contributions only from 
wave vectors in the momentum shell. 

Derivatives of basis vectors can be expanded in the 
local bases 



d,,^ = -1 



(22a) 
(22b) 



in terms of potentials A and B. The arbitrariness of the 
choice of the vectors is reflected in a gauge invariance 
of the potentials i-the gauge transformations are local ro- 
tations around sEJ One can exploit this gauge symmetry 
to show that the potential A^a^i-esponds to higher order 
derivatives of the spin field. tflEil We ignore such contri- 
butions and therefore omit the potential A from now on. 

For the further calculations, x can be considered as 
small since we consider an infinitesimal momentum shell. 
In addition, we treat temperature and disorder, which 
drive the fluctuations of x, as small. In order to expand 
the Hamiltonian in x, we first consider 



-x^BZ^^ + da:^ 

-Xy^X':^ - IxtxtBfa^a + O(x') (23) 



and, ordered by powers of x, 





= vr^+^f)+^r+vp-hO(x^; 


),(24a) 




= j2Br,BT, = j2{d.^r, 

a a 


(24b) 




= j2^d.x:Bra, 

a 


(24c) 




a 


(24d) 


,(2") 


= E {x^xtB^BS, - X^XbBZaB^A 


• (24e) 



An expansion of the energy density for small x gives 



+ Wii'f)[4'^? + 0{x')- (25) 



Substituting this expression back into (pO|), we rewrite 

n = ■Wfroc+Hint (26a) 



^frcc 



Hi, 



(26b) 



d^.E 



-Un4"^)[4M+o{x'). 



(26c) 



Hereby we have separated the "free" and "interaction" 
contributions in a way such that Hint vanishes for x = 
and that Tiftoc contains the bilinear selfinteraction of x- 
Due to the energy contribution Hint the fluctuations 
cannot be integrated out exactly. In analogy to the treat- 
ment of the pure system, we apply standard perturbation 
theory to Hint- lu principle, this can be done in a system- 
atic way at low temperatures, where x ^ T^/^. Aiming 
at the analysis of the stability of magnetic order for low 
temperature and weak disorder, we retain only the renor- 
malization effects to first order in Timt- Integration over 
X leads to an infinitesimal renormalization of the Hamil- 
tonian 



dH — (Hint) 



O(Hfnt) 



(27) 



where the average is over the fluctuations of x weighted 
by Hfroc only. 



Using the averages 

{^^x'^^,x!) 



KfA^ 



one finds 



0, 

0{n) 



{4P = -iN-2) 



(28a) 
(28b) 



(29a) 
(29b) 

(29c) 
(29d) 



Separating the flow of and of by the requirement 
R[{0) = as in Eqs. ( pj| ) for the bare quantities, we 
finally obtain the flow equations 



{d - 2)Kt 



N ~2 
~Kpd 

N - 



^i?^'(0), (30a) 
Kp 



-i^.(^)^.i^.(^)-(2.^ 



(30b) 
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Terms explicitly proportional to the number of replicas 
n have been dropped. A rescaling r — > e'^^r of lengths 
has been included in order to keep the value of the cutoff 
fixed. 

For arbitrary temperature, the flow of K can be inter- 
preted as a flow of the spin stiffness at fixed tempera- 
ture. Using the dimensionless stiffness ji := Ki/K, the 
dimensionless temperature t := T/J = 1/K, the rescaled 
field (j) :— AKip (we recall that K and A are unrenor- 
malized quantities), and Ri{(})) :— Re{'il')j we rewrite the 
flow equations as 

rl N - 2 4A^ ' 

{d-2)n - ——t - _i?^'(0),(31a) 



d£ 



2A 
+ — 



(2 



N ■ 



31 



2 t 



(31b) 



In this form, temperature and disorder strength A appear 
as explicit parameters. The flow equations (|3^) have to 
be solved with the initial conditions 



Jo 



1-pA, 



i?o(0) = In 



1 — p + p e 



0/2 



(32a) 
(32b) 



Corresponding to the neglect of higher orders of Timt 
in Eq. (p7|), these flow equations contain renormalization 
effects only to the leading order in temperature and dis- 
order. Higher orders of perturbation theory would cer- 
tainly generate higher order contributions as well as a 
more complicated functional form of the Hamiltonian. 

Anticipating that i?"(0) > is preserved under the 
flow, Eq. (31a) shows that both thermal fluctuations 
and disorder tend to reduce the effective stiffness. 

In the given order, several features of the flow equa- 
tions are remarkable. The initial function Ra{(j)) depends 
only on p (cf. Fig The crossover from the linear 
regime at small (j) to the linear regime at large </> oc- 
curs at 0* — 21ni^, where the curvature Rq{4>) = 

l/16cosh^[(0 - (/)*)/4] has its maximum R'^^ct)*) = j^. 
The flow equation of the stiffness couples to disorder only 
through -Ro(O). For A = 0, this coupling vanishes, while 
Ri{(j)) = for p = 0. The flow equations depend on N 
only at finite temperatures (unlike the flow equations of 
Ref. ^4|). Thus, the properties at zero temperature are 
expected to be independent of N . However, one has to 
keep in mind that the equations apply only to spins with 
a continuous rotation symmetry {N > 2) and that the 
renormalization scheme ignores the effects of topological 
defects, which are known to be particularly important for 
N = 2. 



IV. RESULTS AND DISCUSSION 

In order to determine the large-scale properties of the 
model, we numerically integrate the coupled flow equa- 



15.0 



10.0 



o 

(DC 




40.0 



FIG. 1: Plot of the initial function Ro{(l>) for the values 
p = 10~^, 10"^, 10"", and 10"^. Note that this nonanalytic 
function is very small up to 0* = 21n where it assumes 
slope that is approximately independent of small p. 



tions (|T]) for ji and ReX't') from ^ = to large i. If 
je converges for £ — > cxd to a finite value much larger 
than t, the system is in an ordered state with finite 
renormalized stiffness and infinite correlation length, i.e. 

(S'(r)S'(r')) decays slower than exponentially with the 
distance |r — r'|. If, on the other hand, ji becomes of 
order i on a finite scale £ = £*, the system is disordered 
and we identify the correlation length as 



(33) 



with some constant A of the oider of the lattice spacing. 
Following previous references, Q'O we specifically define 
this scale from 



t 

2^' 



(34) 



The critical disorder strength Afm separating a ferro- 
magnetically ordered phase from a disordered phase can 
be identified with the line where ^ diverges. Unfortu- 
nately, this criterion does not allow for a distinction be- 
tween true long-range order and quasi-long-range order 
in the ferromagnetic phase. 



A. Zero temperature 

At T = 0, we find the transition line Afm(p) as shown 
in Fig. 1^. This transition has the following features: 

Regime of small p. For dilute disorder, p ^ 0, we find 
a slow but unbounded increase of Afm(p)- In the range 
10^^ < p < 0.1 the line follows roughly the relation [cf. 
the dashed line in the inset of Fig. B 



Afm(p) « 0.606 -h 0.648 In- 



(35) 
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FIG. 2: Two-dimensional phase diagram at T = in a lin- 
ear and semi- logarithmic plot. The solid line represents the 
transition line Afm(p) from the magnetically ordered to the 
disordered phase. In the inset the dashed line represents the 
linear fit (BSl) and the three filled squares visualize the p oints 
with A = 3.7 with p = 0.01, 0.015 and 0.02 used in Sec. |IVC 
for comparison with experiments. 



FIG. 3: Plot of £,{x) for r = with A = 3.7 and A = 22A 
in order to fit the experimental data (filled squares with error 
bars) of Ref. The inset shows a double-logarithmic plot of 
^ as a function of x — Xc for A = 3.7 with Xc = 0.019062. Open 
circles connected by a line represent ^ calculated numerically 
from the flow equations. The dashed line is the best linear fit 
with slope = -0.50007. 



However, in the limit p — )■ 0, Afm(p) appears to increase 
slower than logarithmically. 

This finding imphes order below a finite defect con- 
centration even for arbitrarily strong defect bonds. This 
behavior is in disagreement with Refs. |l^ and |4[ where 
disorder (with A > Agingio) was suggested to destroy or- 
der for infinitesimally small p. As already stated in Sec. 
II D, we believe that the effects of disorder are overes- 



timated in both previous approaches because of special 
assumptions about the nature of the bivectorial field fi. 

For the special case N = 2, where this disagreement 
persists, further references can be included in the com- 
parison. (We assume here that XY and Heisenberg sys- 
tems should have a similar phase diagram at T = since 
our lowest-order flow equations are independent of N. 
Topological defects are ignored in our work as well as in 
Ref. Their presence may further reduce ^.VJ'he way 
how disorder is treated by Cherepanov et aLl3 implies 
that their analysis actually describes an XY model with 
random phase shifts. Various recent work (see e.g. p^ ) 
has provided evidence that quasi-long-range order should 
exist for weak disorder (even in the presence of vortices). 
This observation is consistent with our flow equations b 
contradicting the flow equations of CherepanoV|-ef al. 



From their numerical data, Gawiec and Grempelcrl argue 
for a disordered phase for A > Agingic and p > 0, i.e., 
for Afm(O) = 2. While this conclusion is again in dis- 
agreement with our ijesjult, the numerical data are not: 
Gawiec and Grempelc3 present data for A = 4 as the 
only value with A > 2 and they demonstrate the absence 
of order only for p > 0.02. For A = 4, we flnd order at 
very small p < 0.0062, which actually is not excluded by 
the numerical data (cf. Fig. 17 of Ref. |2^). 



Regime of large p. For p ^ 1~ we find that Afm(p) 
1+ in a smooth way. While this is qualitatively correct, 
a horizontal segment with Afm(p) = 1 beyond the per- 
colation transition at p = ^ is absent. However, in this 
range we find a Afm(p) which is only slightly larger, i.e., 
this is a quantitative effect which may be attributed to 
the continuum approximation as argued below Eq. ([2^). 
In addition, the regime near A = 1 and p ^ ^ is a regime 
of strong disorder in the sense of Sec. II C since >1. 



There our flow equations are not quantitatively reliable 
also because of the lowest-order truncation. On the other 
hand, the transition at small p should be well described 
since Eq. (||) is satisfied, pAp^{p) for p — > 0. 

When the transition line is approached at finite p, the 
correlation length displays a divergence 



e(A,p)~[A-AFM(p)]^ 



(36) 



From the numerical integration of our flow equations we 
find with a mean-field like exponent i/ ~ 0.500(1) [cf. the 
inset of Fig. ||. 



B. Finite temperatures 



Temperature enters the fiow equations in two places 



where it could lead to contrary effects. In Eq. (31a), an 



increase of temperature leads to a faster renormalization 
of the spin stiffness to smaller values [ignoring the tem- 
perature dependence of i?"(0)]. On the other hand, in 



Eq. (31b) temperature tends to suppress i?, which might 



in turn reduce the efficiency of disorder in suppressing 
the spin stiffness. However, from our fiow equations, we 
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always find that thermal fluctuations reduce stiffness and 
therefore also ^. 

More precisely, ^ decreases monotonously with increas- 
ing temperature. ThuSy-a reentrant temperature depen- 
dence as found by GIEj is absent in the present treat- 
ment. We attribute this discrepancy to the fact that GI 
keep the length of the bivector fi fixed, whereas its typi- 
cal length should vanish at low T according to Eq. (12). 
Thereby, with decreasing temperature, the strength of 
disorder is increasingly overestimated, giving way to an 
apparent reentrance. 



C. Comparison with experiments 

We now turn to check the consistency of our theory 
with measurements on La2-£!;Sra;Cu04. To allow for a 
comparison of our results with those of Cherepanov et 
al.ci rWie refer to the same experimental data by Keimer 
et alEB for x = 0.02, 0.03 and 0.04. For the moment we 
assume that interlayer couplings and spin anisotropies 
can be neglected and come back to this issue later on. 

The comparison of our results for ^ with experiments 
involves four parameters, p, A, J, and A. As already 
stated above, our model parameter p = x/2 is directly 
related to the dopant concentration x. ^ is a length scale 
of the order of the lattice spacing. In principle, this pa- 
rameter can depend on temperature and on disorder itself 
(compare the discussion in Rcf. |l^). Such dependences 
could modify the function ^(x,T) in a subdominant way 
and would involve additional assumptions and parame- 
ters. We refrain from including such dependences for the 
purpose of the subsequent semiquantitative comparison. 

At T = 0, J does not enter the flow equations since 
it simply sets the energy scale and A is the only un- 
known model parameter which enters the numerical cal- 
culation of £* . From the consideration of the sujjierex- 
change across a defect bond one expects A » i.e., 
a value clearly above the canting instability. The flnite- 
ness of the measured values of ^ for a; > 0.02 implies 
that the data points lie in the disordered phase, i.e., 
A > Afm(p = 0.01) « 3.67. 

We have determined values of the parameters A, A, 
and J from the requirement that our theoretical values 
for ^ should be consistent with experimental data. In 
view of the given experimental errors and the approxi- 
mate nature of our RG calculation, we found a satisfac- 
tory agreement for 



0.040 



A = 3.7, 

J = 240 K, 
A = 5.8a, 



(37a) 
(37b) 
(37c) 



where a — 3. 8 A is the lattice spacing. In Fig. ^ we 
compare our theory with data for ^ as a function of x 
at T = 0. The inset shows shows for A = 3.7 a double- 
logarithmic plot of ^(x) which reveals the mean-field like 
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FIG. 4: Plot of CHT) for x = 0, 0.02, 0.03, and 0.04. 
Symbols with error bars are experimental data from Ref. 
Open circles connected by lines represent our theory. 



divergence of ^ according to Eq. (^ near the order- 
disorder transition. 

The temperature dependence of ^ is compared in Fig. 
^. Keeping in mind the strong fluctuations of the exper- 
imental data, they can be considered as consistent with 
our analysis. However, the theoretical dependence of ^ 
on T andpO does not quantitatively confirm the empirical 
formulaa'ES 



(38) 



which would imply that in Fig. ^ the curves for different 
p should differ only by a vertical shift. In contrast, we 
find that thermal fluctuations lead to a stronger increase 
of for smaller p. 



D. Coupled layers 

We now address the effects of a very weak interlayer 
coupling J^. For La2Cu04, a ratio J-^/J = 5 x 10~^ 
was detepmined from experiments. eI From a simple scaling 
analysis,Ej one immediately obtains the flow equation 



(39) 



which shows the strong relevance of this coupling. From 
the condition that the interlayer coupling becomes com- 
parable to the intralayer coupling, = J, one can fix a 
dimensional crossover scale 



(40) 



2 J-L 



Comparing this scale with the correlation length (|3^ 
obtained in the absence of the layer coupling, one ex- 
pects that the coupling actually is irrelevant for i* < £ 
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FIG. 5: Plot of our numerical result for Tj^{xi-{open circles 
connected fev a line). For comparison, squaresJSJ trianglesfj 
and crossea^j represent experimental data (cf. the main text). 

where 2D fluctuations on small scales renormalize to 
zero. On the other hand, for £* > £^ fluctuations be- 
come three-dimensional on scales larger than £^ . Per 
definition of £*, the exchange coupling on the scale is 
still large compared to temperature such that magnetic 
long-range order should be stable. Therefore, the loca- 
tion of three-dimensional ordering in the (x, T) plane can 
be determined from the implicit condition 

r(x,T) = 1 In (41) 

We have evaluated this condition numerically for / J — 
5 X 10~^ and the parameter set (p7|). The resulting Neel 
temperature is plotted in Fig. |^ as a function of doping. 
The transition temperature is normalized by its value 
Tn(0) « 300K in the absence of disorder, which essen- 
tially reflects the value of J. While the precise value of J 
may vary to some extent with the employed fltting proce- 
dure for the parameter set, the shape of normalized tran- 
sition line is very robust. At T = 0, the critical disorder 
strength x^d is increased by the interlayer coupling only 
slightly over X2Ti ~ 0.1906, (xsd - X2d)/x2d ~ 10~^, be- 
cause of the extremely slow divergence of ^ for x a; . 

In comparison to the theory in Ref. |lj, we find that 
T^{x) decays less abruptly near x^u- rWe oppose our re- 
sults to data from Hall measurements^ (squares in Fig. 
H), susceptibilitAr measurementsEj (triangles), and PAC 
measurementa23 (crosses). The overall agreement is sat- 
isfying, although the experimental data partially suggest 
that disorder is less effective in reducing the transition 
temperature. This tendency may be attributed to the 
fact that the samples were partially oxygen-doped (at 
higher temperature, oxygen is not quenched and less ef- 
fective in generating spin frustration) as well as the fact 
that we have neglected the easy-plane spin anisotropy, 
which also tends to stabilize the magnetic order. 



V. CONCLUSIONS 

In this article, we have reexamined a classical model 
for A^-component spins (A^ > 2) with random exchange 
couplings. Thereby we have chosen an approach comple^ 
mentary to previous studies by Glazman and loselevichllS 
and by Cherepanov et aZo Special care has been taken to 
preserve the quenched nature of disorder and the global 
spin rotation symmetry. 

Our analysis involves approximations which we briefly 
summarize, (i) The originally discrete spin system is rep- 
resented in a continuum formulation. Thus, features re- 
lated to the speciflc lattice structure - such as the lo- 
cation of the percolation transition for A = 1 - cannot 
be captured quantitatively, (ii) As in the pure case, the 
renormalization scheme accounts only for interactions be- 
tween spin waves, ignoring the role of topological defects. 
Thus, the degree of magnetic order may be overestimated 
(recall that for A^ = 2 the pure system erroneously ap- 
pears to be ordered at all temperatures if vortices are 
neglected) . (iii) The flow equations are truncated to low- 
est order in temperature and disorder. In principle, the 
analysis could be extended to higher orders. In practice, 
this extension is hampered by a much more complicated 
functional form of the Hamiltonian which is generated 
during the flow. 

Due to the nonanalytic nature of the cumulant function 
R, we found it necessary to develop a functional renor- 
malization group, i.e., to keep track of arbitrarily high 
cumulants of the disorder distribution. To the best of 
our knowledge, functional flow equations for disordered 
spin systems have been considered previously only for dif- 
ferent types of disorderp-jpi papticular for random fields 
and random anisotropy.EZlL3'P^tr°l 

The flow equations (|T]) are the central result of the 
analytic part of this work. For comparison, the analysis 
of Cherepanov et aZ.,Lj which employs approximations 
corresponding to the ones listed above, is restricted to a 
single disorder parameter that corresponds to the lowest 
cumulant of the disorder distribution. 

As main physical result from our RG analysis, we find 
that the two-dimensional spin system can be magneti- 
cally ordered at T = in the presence of a sufficiently 
small but finite concentration of arbitrarily strong defects 
(A > Asingie)- At this place it is worthwhile to recall that 
we have identified magnetic order from the length scale 
where the spin stiffness is renormalized down to the scale 
of thermal fluctuations. In the absence of explicit calcula- 
tions for the spin-spin correlation function it is natural to 
assume that this scale coincides with the magnetic corre- 
lation length ^. Such a calculation would be desirable in 
order to clarify whether the ordered phase (with ^ = oo) 
has quasi-long-range order or true long-range order. We 
have determined ^ by a numerical integration of our flow 
equations for the strictly two-dimensional system as well 
as for weakly coupled layers. In the flrst case our zero- 
temperaturej)hase diagram is consistent with numerical 
simulations .E3 In the second case, the calculated depen- 
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dence of ^ on temperature and on disorder strength is 
in good agreement with measurements on cuprates. In 
both cases the comparison was restricted to finite length 
scales given by the computationally manageable system 
sizes or by the experimental error bars, respectively. 

Nevertheless, concerning the question of whether mag- 
netic order is stable against bond disorder in two dimen- 
sions on largest scales, our positive answer, shared, by 
Ref. 15, disagrees with previous negative ones.EJtJ In 
the end of Sec. 11 D we have given specific reasons why in 
our opinion the previous answers cannot be considered as 
final. In view of the complementary approaches and ap- 
proximations involved, this question has to be considered 
as an open one that calls for additional future research. 



for all temperatures. 



APPENDIX B: CONTINUUM LIMIT 

We show that spin frustration mechanism is well cap- 
tured in the continuum representation of our model. To 
this end, we rederive the canting instability threshold 



srlc 



As 

bond located at r 
Hamiltonian reads 



d for a single defect bond. For a single defect 
and oriented in direction e,- the 



H = 



J 



(Bl) 
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APPENDIX A: DUALITY 

For the bimodal distribution (||) with A > 1, our 
system consists of a fraction p of "defective" antiferro- 
magnetic bonds and a fraction 1 — p of "regular" fer- 
romagnetic bonds. The relative strength of bonds is 
\Jafm/ Jfm\ = 1 — a. From a dual point of view, one may 
say that the system consists of a fraction p := 1 — p of 
"defective" ferromagnetic bonds and a fraction 1 —p = p 
of "regular" antiferromagnetic bonds. Since for classi- 
cal spins thermodynamic properties are invariant under 
flipping one sublattice and reversing the sign of the ex- 
change coupUng ( Jafm =: -Jfm and Jfm =: - Jafm),^ 
the system is equivalent to a system with a fraction 
p := 1 — p of "defective" antiferromagnetic bonds and 
a fraction 1 —p = p oi "regular" ferromagnetic bonds of 
relative strength |Jafm/^fm| = 1/(1 ^ A) =: 1 — A.E£l 
Thus duality provides a mapping 



J 
P 
A 



J = (1 - A) J, 
p = 1 - p, 
A 



(Ala) 
(Alb) 



A 



We introduce a canting field x(i")i X^(i") ^ 1j perpendic- 
ular to a coUinear ground state <So(r) = 5*0 of the pure 
system via 



(B2) 



To show the canting instability we insert (|B2| ) into the 
Hamiltonian and minimize the energy with respect to the 
canting field x after expanding the Hamiltonian up to 
quadratic order in this field. In Fourier space the saddle- 
point equation reads 



^1 " ^% / ''i^^- 



(B3) 



Multiplying equation ( p33| ) with qi and integrating over 
q we get 



with 



. _ A 



fit := / ^iXk- 

/k 



(B4) 



(B5) 



This self-consistency condition on pii implies 

A = Asingic = d (B6) 



for a non- vanishing solution of the saddle-point equation 



1' 
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